Electronic states and Landau levels in graphene stacks 
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We analyze, within a minimal model that allows analytical calculations, the electronic structure and Landau 
levels of graphene multi-layers with different stacking orders. We find, among other results, that electrostatic 
effects can induce a strongly divergent density of states in bi- and tri-layers, reminiscent of one-dimensional 
systems. The density of states at the surface of semi-infinite stacks, on the other hand, may vanish at low 
energies, or show a band of surface states, depending on the stacking order. 
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I. INTRODUCTION 



Recent experiments show that single layer graphene, 
and stacks of graphene layers, have unusual electronic 
propertiesii^, that may be useful in the design of new elec- 
tronic devices-''*'^. Among them, we can mention the Dirac- 
like dispersion relation of a single graphene layer*"-^, the chiral 
parabolic bands in bilayers that lead to a new type of Quantum 
Hall effect^, or the possibility of confining charge to the sur- 
face in systems with a few graphene layers''. While the study 
of graphene multilayers is still in its infancy'", the current ex- 
perimental techniques allow for the extraction/production of 
multi-layers with the accuracy of a single atomic layer. The 
ability of creating stacks of graphene layers can provide an 
extra dimension to be explored in terms of electronic prop- 
erties with functionalities that cannot be obtained with other 
materials. 

The nature of the stacking order in graphene multilayers 
has its origins in the single graphene plane, a two-dimensional 
(2D) honeycomb lattice with two inequivalent sub-lattices, A 
and B (two atoms per 2D unit cell). The staggered stacking 
occurs in highly oriented pyrolytic graphite (HOPG) where 
the graphene layers are arranged so that there two graphene 
layers per three-dimensional (3D) unit cell in a ABAB ■ ■ ■ 
sequence (see Fig. Q. In this case, one of the atoms in one 
of the planes (say B) is exactly in the center of the hexagon 
in the other plane. Hence, only one of the sublattices (A) is 
on the top of the same sublattice in the other layer Neverthe- 
less, this is not the only possible stacking observed in these 
systems". Rhombohedral graphite with stacking sequence 
ABC ABC ■ ■ ■ has three layers per 3D unit cell (in the third 
layer the sublattice B becomes on top of sublattice A of the 
second layer). Furthermore, stacking defects have been re- 
peatedly observed in natural graphitic samples'^'^, and also 
in epitaxially grown graphene filmsli. The richness of possi- 
ble stacking configurations is due to the weak van der Waals 
forces that keep the layers together Furthermore, the elec- 
tronic structure of the conduction band in graphene and bulk 
graphite is well described by a tight binding model which in- 
cludes hopping between the tt orbitals in Carbon atoms, ne- 
glecting the remaining atomic orbitals which give rise to the 
a bands in g]-aphit©ii*i^'i^ii2iiS*i2i22. It is known that the low 
energy electronic structure depends on the stacking order in 



bulk s amplest 

In the following, we analyze the electronic structure and 
Landau levels of graphitic structure with different stacking 
orders. We use a minimal tight binding model that only in- 
cludes interlayer hopping between nearest neighbor Carbon 
atoms in contiguous layers, and the k • p expansion for the 
dependence on the momentum parallel to the layers. Within 
the k • p approximation, the two inequivalent corners of the 
Brillouin zone can be studied separately, and we will analyze 
only one of them. We do not consider the effects of disorder, 
already discussed in other publications'"-^^, but focus instead 
on the analytical solution of the model that provide informa- 
tion on the unique properties of these systems. 

The model used here does not take into account hopping 
which may be relevant in order to describe the fine details of 
the band structure, such as the trigonal distortions, that are 
also difficult to estimate using more demanding local density 
functional (DFT) methods. Moreover, recent angle resolved 
photoemission experiments in crystalline graphite^** show that 
those effects are very small. On the other hand, the methods 
used here, based on a mapping onto a set of simple nearest- 
neighbor one-dimensional (ID) tight binding Hamiltonians is 
quite feasible, allowing us to study a variety of situations, with 
and without an applied magnetic field. It is worth noting that 
when the energy scales associated to the electron-electron in- 
teraction and/or to in plane disorder are larger than the hop- 
ping neglected here, the description presented below will be 
a reasonable ansatz for the calculation of the effects of the 
electron-electron interaction. 

The model used is described in the Section HII as well as a 
simple scheme that allows the mapping of the problem of an 
arbitrary number of coupled graphene layers onto a ID tight 
binding model with nearest neighbor hopping only. Section 
Unidescribes the main properties of few layer systems, also in 
the presence of a magnetic field. In Section lTvl we discuss the 
bulk and surface electronic structure of semi-infinite stacks 
of graphene layers, with different stacking order. Section |V] 
contains the main conclusions of our work. 



II. THE MODEL. 

We consider the staggered stacking, where the layer se- 
quence can be written as 1212 • • • , and rhombohedral stack- 




FIG. 1: Top view of the arrangement of the Carbon atoms in two 
neighboring graphene planes in the staggered stacking. 
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FIG. 2: Schematic view of the mapping of the problem of a stack of 
graphene layers onto a set of one dimensional problems. The hop- 
ping between sites are t± and fce'*"*, where (j) — arctan{ky/kx). 
Left: staggered stacking. Right: Rhombohedral stacking. Note that 
the phases (j> can be gauged away in both cases. 



ing, 123123 • • • and combinations of the two. We do not con- 
sider hexagonal stacking, where all atoms in a given plane 
are on top of atoms in the neighboring layers. 111 • • • . In all 
stacking considered, the hopping between a pair neighboring 
layers takes place through half of the atoms in each plane, as 
schematically shown in Fig. [2 

We describe the in-plane electronic properties of graphene 
using a low energy, long wavelength, expansion around the 
K and K' corners of the 2D Brillouin Zone. This description 
uses as only input parameter the Fermi velocity, vp — 3ia/2, 
where i (« 3 eV) is the hopping between tt orbitals located 
at nearest neighbor atoms, and a = 1.42 A is the distance be- 
tween Carbon atoms. We assume that hopping between planes 
takes place only between atoms which are on top of each other 
in neighboring layers. We denote the hopping integral as t± 
(« O.li). The whole model is defined by the parameters vp 
and t±. For instance, for the graphene bilayer (the 3D unit 
cell of the staggered stacking) the Hamiltonian has the form: 
H = ^i^*t(k)iyo(k)*(k), where k = ik^,ky) is the 2D 
momentum measured relative to the K (K') point (we use 
units such that c = 1 = h). 
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(1) 



*^(k) = (cl_i,k 4,i.k cl,2,k 4,2,k) is the electron spinor 
creation operator, and 0(k) = ta.'n~^ {ky / k^) is the 2D angle 
in momentum space. 

In the case of the staggered stacking with 2N planes, the 
full Hamiltonian consists of 2 A^ x 2 A^ matrix with N/2 blocks 
with size 4x4 given by 0. The spinor operator is given by 
creation and annihilation operators for electrons in different 
planes: (y^^'k. where a = A^B labels the sublattices in each 



plane, and n — !,■ 
Green's functions: 



, N labels the planes. We define the 



GS,Mk,i) = -^(^c„,„,kW4,„,k(0)) 



(2) 



where T is the time ordering operator, and its Fourier trans- 
form, G^ ^(k, ijj), that can be used to calculate the properties 
of these systems. 



It is easy to see from Q that G^(|k|, w) {a - 
G^(|k|, w) (a — (3 — B), satisfy the equations: 
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The phase factors e*'^'', e~^'^^ can be gauged away by a gauge 
transformation, leading to a set of equations where the odd 
and even numbered layers are indistinguishable. We can inte- 
grate out the B sites in eq.(|3j, and write an effective equation 

forGA(|k|,c^): 
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which is formally equivalent to the tight binding equations 
which describe a ID chain with hopping energy t± and ef- 
fective on site energy (wF|k|)^/a;. The Green's function for 
the B atoms can be obtained from Gyi(|k|,a;) using the ex- 
pression: 
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A sketch of the resulting ID tight binding model is given in 

Fig.ia 

In the case of the rhombohedral stacking (123123 ■■■'), a 
suitable gauge transformation can be used in order to make 
all layers equivalent, although the two sublattices A and B 
within each layer remain different, leading to the following 
set of equations: 
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(7) 



that determine the properties of the rhombohedral stacking. 

In a finite magnetic field B the above equations can be mod- 
ified with the use of a Peierls substitution k ^ k + e A, where 
A is the vector potential (V x A = B). In this case, the 
system breaks into Landau levels that, for a single graphene 
layer, have energy^^: 



E{s,j) = SUJcy/j + 1, 



(8) 



where j — 0, 1, • • • labels the Landau levels, s = ±1 labels 
the electron and hole bands. 



= \/2i;f/^b = VFV2eB, 



(9) 



is the cyclotron frequency, and Ib is the cyclotron length. 

In the presence of magnetic field we can replace the mo- 
mentum label k by the Landau label j (|k| -^ Vj/^b) and the 
extension of the equations Q become: 
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Note that, in this case, the two inequivalent layers in the unit 
cell cannot be made equivalent by a gauge transformation as 
in(|3}. 

In rhombohedral case, the presence of a magnetic field, (0 
become: 

uG\{^) = t^Gl-\u) + ^VjGliio) , 



UjGUlo) = t^Gl+\L0) + ^^/JTlG\{L0). (11) 

III. STACKS OF A FEW GRAPHENE LAYERS. 

A. Generic energy bands. 

Using the ID tight binding mapping discussed in the pre- 
ceding section, we can write an implicit equation for the 
eigenenergies of a system with N layers with the staggered 
stacking as: 
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where n — 1,- ■ ■ N,so that the dispersion becomes: 



(12) 
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FIG. 3: The Fermi ring: the Fermi sea of a biased bilayer. 



B. Bilayer. 

A number of properties of the present model applied to a 
graphene bilayer can be found in ref. 112^ . We extend those 
results here to the case where there is an electrostatic potential 
which makes the two layers inequivalent. The Hamiltonian 
reads: 
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(14) 
where the difference in the electrostatic potentials in the two 
layers is 2A. For UF|k| ^ A ^ t±, the bands at low energy 
are given by ±ek where: 
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Hence, the bands show an unusual "mexican hat" dispersion, 
with extrema at fcg ^ ^/^f with energy e^H = A — A^ /t']_. 
Near these special points, the electronic density of states di- 
verges as: 
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which is the divergence seen in ID systems. In order to un- 
derstand the reason for this ID behavior consider the situa- 
tion where the chemical potential is above (or below) £„//. In 
this case the Fermi surface looks like a ring (a Fermi ring), 
as shown in Fig. |3l At large radius fco the ring approaches a 
ID dispersion with nested Fermi surfaces, characteristic of ID 
systems. 

In the presence of a magnetic field, the Hamiltonian be- 
comes: 
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where j is the index of the Landau level. For values of j such 
that (wp V7)/^B < A < ij^, we obtain: 
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C. Trllayer. 

1. HOPG stacking (ABA). 

The model used here leads to six energy bands in a trilayer. 
In the absence of electrostatic potentials that break the equiv- 
alence of the three layers, we can use ea. (ll3t . to obtain; 
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There is a band with Dirac-like linear dispersion, and two 
more bands which disperse quadratically near e = 0, as in 
a bilayer The Landau levels, at low energies are: 
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(20) 



Hence, the spectrum can be viewed as a superposition of the 
corresponding spectrum for the Dirac equation, and that ob- 
tained for a bilayer. 

The calculation of the energy levels in the presence of elec- 
trostatic fields which break the symmetry between the three 
layers is more complex, and cannot be performed fully ana- 
lytically. We consider the case when the upper layer is at po- 
tential A, the lower layer is at potential —A, and the middle 
layer is at zero potential. The Hamiltonian is: 
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(21) 



In the limit up |k| ^ A ^ <j^, one can use perturbation theory 
to obtain: 
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This equation, although approximate, describes correctly the 
existence of states at zero energy when wp|k| — A. The dis- 
persion shows extrema at momenta fco oc A/vp with energy 
e^H oc A'^ /t±, leading to a divergent density of states: 
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This expression, as in the case of eq. ( I16> . is typical of ID sys- 
tems. The low energy bands of a trilayer with ABA stacking 
are shown in the upper panel of Fig.l4l. 



2. Rhombohedral stacking (ABC). 

The low energy bands of a trilayer with the ABC stacking 
differ significantly from those for the ABA stacking discussed 
previously. The 6x6 hamiltonian is: 
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Assuming that A -^ t±, only two out of the six n bands lie at 
energies |e| ^ t±. These bands are derived from the orbitals 



at the B sublattice in the two outermost layers, given by the 
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FIG. 5: Density of states for the staggered stacking. Continuous line: 
IinGs(Lj). Dashedline: YroGAiyi)- 



FIG. 4: Low energy bands of a trilayer with an electrostatic field 
which lireaks the equivalence of the three layers (see text). Top: 
HOPG stacking, ABA. Bottom: rhombohedral stacking, ABC. The 
parameters used are vp = l,t± — 0.1, A = 0.04. 



We can obtain the local density of states by integrating 
Gyi(|k|,Ci;) and GB(|k|,a;) over the parallel momentum k. 
This integral can be performed analytically, and we obtain, 
for the bulk; 



first and sixth entry in ea.i24\. One can define an effective 
2x2 hamiltonian: 
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The dispersion is cubic at low momenta, ek ~ (tiF|k|)'^/t^. 

The density of states is D{e) oc t/ /(upe^/"^). The low energ 
bands of a trilayer with ABC stacking is shown in the lower 
panel of Fig. 0. 

The effective 2x2 hamiltonian in the presence of an applied 
magnetic field is, for A = 0: 
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(26) 
There are, in addition, states localized in the outermost lay- 
ers, with Landau index n = 0, and in two layers, with Landau 
index n = 1. It is interesting to note that Landau levels asso- 
ciated to different K points are localized in different layers. 
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The local density of states is given by (see also ref. |lul): 
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The density of states in the bulk of the staggered stacking is 
shown in Fig.|5l 

We can use the equivalence to a ID model for each value 
of k, (|3} and Q, in order to obtain the Green's function at the 
surface layer of a semi-infinite system: 



IV. SEMI-INFINITE STACKS OF GRAPHENE LAYERS. 

In this section we consider the case of a graphene multilayer 
with a surface termination. In the case of the staggered stack- 
ing, we can use eq.Q, and obtain for the Green's function in 
the bulk as: 
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(30) 



for staggered stacking. The local density of states, obtained 
after integrating these expressions over k is shown in Fig.|5] 

The ID tight binding model which describes the electronic 
bands in rhombohedral graphite (stacking order 123123 • • •) 



D((0) 




FIG. 6: Density of states at the surface layer. Continuous line: 
IinGfl(a;). Dashed line: ImGA{^)- 




is given in Q. Using this equation, the Green's function, in- 
tegrated over the perpendicular momentum, k±, is given by: 
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The local density of states can be obtained by integrating this 
expression over k: 
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The material, within this approximation, is a semi-metal, with 
vanishing density of states at the Fermi level at half filling, 
uj = 0. Note that the density of states at low energies is inde- 
pendent of the value of t± . 

In the presence of a magnetic field, we use (lll> . We can 
integrate out the Green's functions for the sites in a given sub- 
lattice and obtain: 
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These equations are formally equivalent to those obtained for 
the wavefunctions of a displaced ID harmonic oscillator: 
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(34) 



with the correspondence: cq ^-> t\/^, ojq <-> i;|,/(£|w), g ^^ 
(i_LWF)/(^Bw). The eigenenergies of (l34l are e,„ == eo — 
g^ jijjQ + mujQ, and we can write the eigenenergies associated 
to ea.lfni as: 
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FIG. 7: Average density of states of a stack of 50 graphene layers 
with rhombohedral order, ABC ABC ■ ■ ■ , embedded into the stag- 
gered stacking, ABAB ■ ■ ■ . The parameters used are: ix = O.leV, 
and B = IT (^b « 25.7nm, vy/1^ ~ 0.025eV). The bands of some 
of the lowest Landau levels of the staggered stacking, calculated for 
the same parameters, are shown for comparison. 



The spectrum, in an applied magnetic field, is discrete, and 
equal to that in graphene, as well as the local density of states 
in the absence of the field, eq. M>2\ . The discreteness of the 
spectrum survives when a stack of rhombohedral graphene is 
embedded into the staggered stacking, as shown in Fig.0 

Using the model described in Section |llj we obtain for the 
projected density of states: 
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(37) 



The real part of these Green's functions has a pole at w = 
for |k| < t±/vF, indicating the presence of a surface state. 
The projected density of states, as function of |k|, is sketched 
in Fig. El 

We can also analyze modifications of the surface, induced 
by electrostatic potentials or local changes of the stacking or- 
der A shift of the topmost layer by a potential e does not 
change qualitatively the projected band structure shown in 
FigEl unless e > t±. On the other hand, a stacking of the type 
1232323 • • • leads to a surface band, which can be shifted by 
an external potential, e, applied to the topmost layer, labeled 1. 
As the structure looks like a perfect staggered stacking beyond 
this layers, and the couplings are local, the Green's function 





FIG. 8: Sketch of the projected density of states at the surface, as 
function of parallel momentum. Left: staggered stacking. Right: 
rhombohedral stacking (the line at e = stands for a dispersionless 
band of surface states). 



in the topmost layer is given by: 
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where the labels A^ B correspond to the two inequivalent sites 
in the topmost layer A sketch of the resulting projected elec- 
tronic structure is shown in Fig. |9] This result is consis- 
tent with the observation of at least two frequencies in the 
Shubnikov-de Haas experiments in graphene multilayers with 
induced carriers at the surface reported in ref. |2J- 



V. CONCLUSIONS 

We have analyzed a minimal model for the electronic struc- 
ture of stacks of graphene planes, which allows us to obtain 
a number of interesting results analytically. In the following, 
we outline some of the most relevant ones: 

• Bilayers and trilayers, in the presence of electrostatic 
fields that break the equivalence of the layers, develop 
van Hove singularities where the density of states be- 
haves in a quasi ID fashion, D{e) oc (e — Euh )^^/^. 

• The Landau levels of a trilayer, in the absence of elec- 
trostatic effects, are given by a set of levels which 
depend on the field in the same way as for a sin- 
gle layer, e(j) = ±(t;F-\/7)/^B, and another set 
which is equivalent to the levels in a bilayer, e(j) = 
±((wp-\/j(j + l))/tj_i'|). A comparison of the two 
sets will allow for an independent measurement of the 
value of t±^ . 

• The surface density of states in staggered stacking van- 
ishes at zero energy at the sites with a nearest neighbor 




FIG. 9: Sketch of the projected density of states of a stack of 
graphene layers with ordering CABAB ■ ■ ■ , in presence of a po- 
tential shift in the topmost layer. The lines outside the continuum of 
states stand for a two bands of of surface states. 



in the next layer This result may explain the differ- 
ence in the images of the atoms at the two sublattices 
observed in STM measurements of graphite surfaces^-. 

• The bulk local density of states of rhombohedral stack- 
ing (ABC ABC ■ ■ ■) vanishes at zero energy, and is in- 
dependent of the value of the interlayer hopping, t^. 
The spectrum of bulk Landau levels is discrete, unlike 
most three dimensional systems, and shows the same 
field dependence as that of a single graphene layer, with 
independence of the value of tj_. Hence, discrete, quasi- 
2D Landau levels may exist in nominally 3D samples, 
helping to explain the experiments in ref. 13011 . 

• The surface of rhombohedral, and of staggered stacking 
with a rhombohedral termination {CABAB ■ ■ ■) has 
surface states, with a well defined dispersion as func- 
tion of the parallel momentum. This result may help to 
explain the observation of Shubnikov-de Haas oscilla- 
tions in graphitic systems with a highly doped surface^. 
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